Dynamics of quantum dissipation systems interacting with bosonic canonical bath: 

Hierarchical equations of motion approach 
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A nonperturbative theory is developed, aiming at an exact and efficient evaluation of a general 
quantum system interacting with arbitrary bath environment at any temperature and in the presence 
of arbitrary time-dependent external fields. An exact hierarchical equations of motion formalism is 
constructed on the basis of calculus-on-path-integral algorithm, via the auxiliary influence generating 
functionals related to the interaction bath correlation functions in a parametrization expansion form. 
The corresponding continued-fraction Green's functions formalism for quantum dissipation is also 
presented. Proposed further is the principle of residue correction, not just for truncating the infinite 
hierarchy, but also for incorporating the small residue dissipation that may arise from the practical 
difference between the true and the parametrized bath correlation functions. The final residue- 
corrected hierarchical equations of motion can therefore be used practically for the evaluation of 
arbitrary dissipative quantum systems. 
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I. INTRODUCTION 

As a fundamental topic in quantum statistical mechan- 
ics, the development of quantum dissipation theory has 
involved scientists from diversified fields of research over 
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The key quantity in quantum dissipation theory is the re- 
duced density operator, p(t) = tr B pr(t), i.e., the partial 
trace of the total density operator over the bath space of 
practically infinite degrees of freedom. The standard ap- 
proach to the exact p(t) is the influence functional path 
integral formalism [l|. However, from both operational 
and numerical points of view, this formally exact formal- 
ism is much limited in comparison with its differential or 
master equation counterpart. 

The attempt to the construction of quantum master 
equation, via the derivative on the exact path integral for- 
malism, was first carried out by Caldeira and Leggett, for 
the quantumdissipation in the high temperature Marko- 
vian limit The calculus-on-path-integral (COPI) 

method has also been used in driven bistable systems, 
together with the so-called non-interacting blip approxi- 
mation or its variations to treat the quantum path corre- 
lations in the reduced dynamics [4|, [8|, [l2j . The first exact 
quantum master equation via the COPI method is con- 
structed for the damped harmonic oscillator systems 31, 
l32l ]. Tanimura and coworkers have recently constructed 
a set of hierarchically coupled quantum Fokker-Planck 
equations for a general Drude dissipation system [2lj . 
without any lon ger the high-temperature approximation 
exploited before [19l. [20j . The deterministic description of 
an exact differential equations of motion (EOM) formal- 
ism can also be constructed via stochastic descriptions 
of quantum dissipation [H, El El El EE El EE EH- 
However, it has by far only been carried out with a cer- 
tain single-mode dissipation model system. The general 
theory of quantum dissipation in terms of deterministic 



EOM remains a great challenge. 

The aim of this work is right at this fundamental issue 
of quantum statistical mechanics. It is to develop an ex- 
act and efficient EOM formalism for a general quantum 
system, interacting with arbitrary bath environment at 
any temperature, and arbitrary time-dependent external 
fields. The only assumption involved is the same as in 
the path integral formalism that the interaction bath sat- 
isfies the Gaussian statistics p], H [l7[. The theoretical 
construction will be made mainly based on the COPI al- 
gorithm developed in Ref. [H A set of so-called influence 
generating functionals will be shown to be crucial in for- 
mulating the exact hierarchical EOM. The present work 
focuses on the canonical bath ensemble case, but can be 
easily extended to the grand canonical case [331 ] . 

The remainder of this paper is organized as follows. 
In SecHH after the description of the general form of 
system-bath coupling Hamiltonian, we reviews the path 
integral influence functional formalism (l7j , together with 
its formal relation to the time-local dissipation superop- 
erator in quantum master equation. In Sec. lIIH we il- 
lustrate with the Drude-Debye model the key technical 
issue on the construction of hierarchical EOM. It is to 
identify a complete set of auxiliary Influence Generating 
Functionals via the COPI algorithm [17|, which will be 
termed as the IGF-COPI for short hereafter. Section 
IIVI turns to general dissipative quantum systems. It in- 
volves a parametrization of general interaction bath cor- 
relation functions that satisfy the fluctuation-dissipation 
theorem; see AppendixfAl for details. Upon identifying 
the complete set of auxiliary IGF-COPI construction for 
the parametrized forms of correlation functions, exact 
hierarchical EOM are followed immediately. The equiv- 
alent formalism in terms of hierarchical and continued- 
fraction Green's functions and memory kernels are given 
in AppendixlBl In Sec. El we further establish the so- 
called principle of residue correction. The principle it- 
self is rooted at the formally exact relation established 
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in Seaim between the Feynman- Vernon influence func- 
tional and the time-local dissipation superoperator. But 
it is now exploited in certain perturbative or nonpertur- 
bative approximate manner to construct an efficient hier- 
archy truncation method. It is also used to incorporate 
the residue contribution, due to the difference between 
the exact and the parametrized bath, into the final the- 
ory. Finally, Sec. IVII concludes this paper. 

II. QUANTUM DISSIPATION IN TERMS OF 
INFLUENCE FUNCTIONALS 

A. Multi-mode system— bath coupling Hamiltonian 

Consider a general quantum system embedded in a 
bath which is assumed as a canonical ensemble in this 
work. The total system-plus-bath composite Hamilto- 
nian can be written in general as 

H T = H{t)-Y,QaF a + h B . (1) 

a 

Here, H(t) and h B denote the uncorrelated system and 
bath Hamiltonians, respectively. The former may be sub- 
ject to a time-dependent coherent field drive. The second 
term in the right-hand-side (rhs) of Eq. ([1]) denotes the 
multi-mode system-bath coupling. It can be generally 
expressed in the multiple dissipative mode decomposi- 
tion form, in which {Q a } and {F a } are the system and 
bath operators, respectively, and assumed to be Hermite 
in this work. 

For the later use, introduce here the Liouvillian C and 
the dissipative coupling mode Q a in the reduced system 
subspace via their actions on an arbitrary operator as 

CO = [H(t),6}, Qad = [Q a ,0}. (2) 

Throughout this work, we set h = 1 and the inverse tem- 
perature (3 = l/(k B T), and denote also dt = d/dt. 

The stochastic canonical bath ensemble average of an 
operator O is denoted as 

(6) B EE tr B (<5/£») = tr B (6e-^)/tr B e-^. (3) 

In the /i B -interaction picture, 

F a (t) ee e^Fae-^, (4) 

for each bath interaction operator in Eq. ([T]) , is assumed 
a Gaussian stochastic process. This is exactly the case 
when the bath consists of a collection of harmonic oscil- 
lators and each F a is a linear combination of the coordi- 
nates and momenta of the harmonic bath oscillators. The 
Gaussian stochastic process is also related to the central 
limit theorem in statistics. As the quantity Q a {F a ) B , if 
it is not zero, can be included in the system Hamiltonian 
in Eq. ([T]), we can set {F a ) B — without loss of general- 
ity. The effects of Gaussian stochastic bath operators are 
then completely described by their correlation functions, 

C ab (t - t) = (F q (*)F 6 (t)} b . (5) 



They satisfy the symmetry and detailed-balance relations 
of the canonical bath 0, UHl > 

C* ab (t) = C ba (-t) = C ab (t - 0). (6) 

We shall be interested in a differential EOM formalism, 
aiming at an efficient and exact evaluation of the reduced 
dynamics, with arbitrary multi-mode non-Mar kovian dis- 
sipation and time-dependent external fields on the sys- 
tem. The key quantity of interest is the reduced density 
operator p(t) , defined together with its associating prop- 
agator U(t, t ) as 

p(t) = tr B p T {t)=U{t,t )p(t ). (7) 

The desired EOM will be constructed (cf. Sec. HV)) via the 
IGF-COPI approach [ljj, starting from the exact path 
integral expression that involves only the initial factor- 
ization ansatz prito) = p{to)p B q . Note that when the 
initial time is set to t,Q — + — oo, this ansatz becomes exact 

am. 

B. Reduced dynamics versus influence functionals 

This subsection summarizes the path integral formal- 
ism of quantum dissipation for the multi-mode system- 
bath interaction [13]. Exploited explicitly is only the 
Gaussian statistical property, the essence of linear har- 
monic bath coupling with arbitrary system operators 
{Qa}- Unlike the EOM formalism that can be expressed 
in the operator level, the path integral expression goes 
with a representation. Let {|a)} be a basis set in the 
system subspace. In the a-representation, Eq. ([7]) reads 
[setting a = (a, a') for abbreviation] 

p(a,t) ee p(a,a',t) = Jda U(a, t;a ,t )p(a ,t ). (8) 

The reduced Liouville-space propagator reads in terms of 
path integral as [l| 

U(a,t;a ,t )= / Va.e lS[a] T[a\e- lS[a ' ] . (9) 

Ja [t ] 

The effect of system-bath interaction on the reduced 
system dynamics is described by the Feynman- Vernon 
influence functional which will be elaborated soon. 
In Eq. ([9]), S[a) is the classical action functional of the 
reduced system, evaluated along the path a{r), with 
the constraints of two ending points a (to) = ao and 
a(t) — a being fixed. In the absence of bath interaction 
(J- = 1), the dynamics would be completely coherent; 
i.e., dtU — —iCU; if T = 1. 

Consider now the key quantity, the bath interaction 
induced influence functional T . Traditionally, its ex- 
pression is derived by adopting a single-mode system- 
bath interaction model P, m which the bath h B is 
assumed to consist of a set of uncoupled harmonic os- 
cillators {qj} and the system-bath interaction assumes 
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the form of H' = —QF — ~QJ2j c j1j' rather than the 
multi-mode decomposition as the second term of Eq. |T]). 

In connection to the later development of EOM formal- 
ism, we denote a = (aa') for a pair of dissipation modes 
hereafter and introduce [cf. the Eqs. (5) and (6) of Rcf. 

m 

Q a (f, {a}) ee Q aa ,(t; {a}) - Q' aa ,(t; {a'}), (10) 
where (noting that Q a = Q aa , and C a = C aa >) 



Q a (t;{a})= / dTC a (t-T)Q a ,[a(T)], (11a) 
)' a (t;{a'})= [ drC* a (t-r)Q al [a'(r)}. (lib) 

J t 



Denote also 



Q a [cx{t)]=Q a [a(t)]-Q a [a'{t)]. 



(12) 



It is in fact the Q a -commutator [cf. Eq. ^ the second 
identity] in the path integral representation, as it depends 
only on the fixed ending points. The final expression of 
the influence functional reads [T3l 



J- [a] ee cxp | — J dr TZ[t; {a}] 



(13a) 



with 



K[t; {a}] ee Qa[a(t)] Q«(t; {a}). (13b) 



The above relations will be used together with Eqs. (fTU)l- 
(fT2")) in the following sections to develop the desired EOM. 



C. Influence functional versus dissipation 
superoperator 



Note that in Eq. (|T3|) the conventional influence phase 
functional is now expressed in terms of its time integrand 
1Z. The latter is in fact the time-local dissipation super- 
operator lZ(t) in the path integral representation, as the 
time derivative of Eq. ([9]) with Eq. (p~3|) leads to 



or equivalently 



d t U = -i£U-K{i)U, 



-iCp — lZ(t)p. 



(14a) 



(14b) 



We can therefore call 1Z of Eq. (|13b[) as the dissipation 
functional. As inferred from Eq. (|13b|) with Eqs. (fTUl) and 
(fT2")l . it may be symbolically expressed in the operator 
level as 



a 



(15) 



Here, Q a denotes the operator form of Eq. (| 1 1 a,|) . How- 
ever, the explicit operator-level expression for 1Z (or Q a ) 
is generally not available. In the presence of a time- 
dependent external field, the only case being of analyt- 
ical expression of 1Z (or Q a ) is, to our best knowledge, 
the driven Brownian oscillator system [161 ]. 

The formal relations, Eqs. (|13p - (fT5"|) , also highlight 
where the difficulty is in the exact evaluation of quan- 
tum dissipation; all relate to the memory-containing Q- 
or equivalent Q-functionals. In a certain sense, the EOM 
formalism to be presented in detail soon (cf. Sec. lIIII 
and Sec.lIVp is a mathematical construction that hier- 
archically resolves the "history" containing in the Q- 
functionals. This issue will become evident in the coming 
sections. Moreover, these formal relations may also be 
directly useful for the residue correction of the final for- 
malism, due to either hierarchy truncation, or the small 
difference between the exact but complicated bath corre- 
lation functions and the parametrized ones (cf. Secly]). 



III. HIERARCHICAL EQUATIONS OF 
MOTION: CORRELATED DEBYE DISSIPATION 

Before presenting the formalism for arbitrary dissipa- 
tive systems, let us consider in this section the simplest 
multi-mode dissipation case, the Drude-Debye model, in 
which 



C a (t>Q) = C aa ,(t)= Va e 



-Ja.t 



(16) 



The parameters j a = j aa > — j a > a are real, while r/ a = 
Vaa' = Va'a are complex, and they are the same as the 
Drude parameters 7™ and 77° in the next section, where 
the general case is studied. Despite this fact, this section 
provides with clarity the basic ingredient of the IGF- 
COPI approach to the desired hierarchical EOM. 

The hierarchy construction starts with the time deriva- 
tive on the propagator U [Eq. @] of the primary inter- 
est. The time derivative on the action functional parts 
contributes to the coherent dynamics of —iCU, and thus 
can be included into the final EOM. We shall show in 
the following that the hierarchy generation stems from 
the time derivative on the evolving influence functionals 
in each tier of the construction. 

Consider first the time derivative on the influence func- 
tional of primary interest [Eq. (fi"3"|) ]. 



dtF[a] = -[Y,Qa[ct(t)]Qa(t;{ct})\f[a 



(17) 



We shall hereafter omit the explicit path integral vari- 
ables dependence whenever it does not cause confusion. 
As results, we recast Eq. (JTTJ) as 



iQc 



T '= -i^QaTa. (18) 



The last identity introduces the auxiliary influence func- 
tionals (AIFs), 



(19) 
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In contrast to T whose leading term is 1, the first tier 
AIFs {Fa} are of the second-order in the system-bath 
coupling as their leading terms. The hierarchy to be 
constructed will go naturally with increasing the order 
of system-bath coupling. 

Consider now the time derivative on the first tier AIFs 
[Eq.®]. 

d t T a = -i{d t Q a )T Qb(-iQa)(-iQb)] T. (20) 



The second term, which arises from the derivative on 
the primary T [cf. Eq. I|18p]. introduces now a set of new 
AIFs, 



T ab = {-iQ a ){-iQ h )T. 



(21) 



The leading terms in these second-tier AIFs are of the 
fourth-order system-bath coupling. 

The first term of Eq. ([20]) involves the time derivative 
on Q a [Eq. (HOI) with Eqs. flTT}]. From Eq. ([TTa|) 



d t Q a = C a {0)Qa'[a(t)] + Qo(t; {a}), 



(22) 



with 



Qa(t;{a}) 



dTC a {t-T)Q a ,[a{T)}. (23) 



Apparently, Q a leads to Eq. (j!?0|) a non-hierarchy term in 
general, unless a specific form of C a (t) such as Eq. (fT6|) is 

considered. In the present case of study, Q a = — 7oQa, 
leading to 

dtQa =iC a -JaQa- (24) 

Here [noting that C o (0) = r] a from Eq. (fl"6"j)] 

C a = -i{VaQa>Ht)] - VaQa>[a'(t)}}, (25) 

which depends only on the fixed ending points. Substi- 
tuting Eqs. jH]) and flU} into Eq. ((2Qj) leads to 



d t T a = C a T - JaF a - i ^ Qbfab 



(26) 



We are now in the position to complete the IGF-COPI 
approach to the hierachical EOM for the multi-mode 
Drude-Debye dissipation of Eq. (fT6)l . First of all, we no- 
tice that due to the mathematical group nature implied 
in Eq. (fTgj) and Eq. (pM]) , (— iQ a ) constitutes the single in- 
fluence generating functional for each pair of the Drude- 
Debye modes. The involving AIFs can be gcncrically 
expressed as [cf. Eqs. (|T5|) and pijl] 



(27) 



Here n = (n a ,n , • ■ • ) consists of a set of nonnegative in- 
tegers. Denote also the index-set, = (n a ± 1, n D , ■ • • ), 



that deviates from n only by changing the specified n a to 
n a ± 1. Note that the total number of nonnegative inte- 
gers in the index set n is the same as that of the non-zero 
system-bath coupling mode pairs. 

The time derivative of T n can be carried out by using 
Eq. (f24|) and the first identity of Eq. (fT8|) , resulting in 

dtT n = -(£n a j a )J r n +J2( naCa:F ^- iQ " :F n +). (28) 

a a 

Define now the auxiliary propagators by [cf. Eq. |9])] 

pa 

U n {a,t;a ,t ) = / Vae lS[a] T n [a]e~ lS[a ' ] , (29) 

which also define the related auxiliary density operators, 
p n (t)=U n (t,t )p(t ). (30) 
Equations ([28]) can then be recast into 

p n = - (iC+^riaja) p n +^(n a C a p n - -iQ a p n +) ■ (31) 



The involving Q a and C a , which in Eqs. (|28|) were given 
by Eq. ((TJJ) and Eq. (|23|) . respectively, at the fixed ending 
points of ck(t) = a in the path- integral representation, 
are now defined in the operator level by Eq. ([2]) and 



C a O = -i(l]aQa'd - T] a OQ a ' 



(32) 



Each of the summations in the rhs of Eqs. (I3ip runs over 
all non-zero coupling mode pairs a = (aa'). The initial 
conditions to Eqs. (|3"Tj) and the methods of infinite hier- 
archy truncation will be discussed later together with the 
general dissipation systems; see comments after Eqs. (|50|) 
and in Sec.lVl 

Equations (I31j) generalize the previous work on the 
sin gle- mode Drude-Debye dissipation [13, [H, H3, HH, HI, 
[23l. |24|. It is noticed that at the second or higher tier, 
which is of the fourth or higher order in the system-bath 
coupling, contributions from different dissipative modes 
are no longer just additive. Some interesting phenomena 
such as co-tunneling [34[ and co-dissipation (e.g., coop- 
erative T2-decoherence and Ti-relaxation processes) will 
therefore be anticipated and investigated elsewhere. 



IV. HIERARCHICAL EQUATIONS OF 
MOTION: GENERAL NON-MARKOVIAN 
DISSIPATION 

A. Non-Markovian bath via parametrization 

It is evident now that the construction of hierarchical 
EOM for a general dissipative system should involve a 
proper parametrization scheme for C aa '(t). It is required 
by the IGF-COPI structure that all involving dtQ a terms 
be contained within the hierarchy; cf. Eqs. (f2"2"| - (|24p and 
comments there. On the other hand, the relations of 
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Eq. ([6]), or more precisely the fluctuation-dissipation the- 
orem (FDT) should also be observed. 

In this work, we adopt a FDT-preserved parametriza- 
tion scheme, in which the bath correlation functions 
C a (t) for a general system at any temperature are com- 
pletely characterized by a set of real parameters; see 
details in Appendix[X] In particular, the parameters, 
{7",7",w"; k = 0,---,K} that will explicitly enter the 
final EOM, together with the Matsubara frequencies 
{7 m = 27rm//3; m—1,- ■ -,M} are all positive, except u>q = 
0. For the latter use, we denote also ui' k a = d k olo + w fc- 

The final FDT-preserved bath correlation functions as- 
sume the following form [cf. Eq. (|A4|) where M — * oo ] . 

2K+1 M 

C a (t>0)=T&e-TZ*+ ^^^W + ^C^™' (33) 

j—Q m—1 

Here, j m — 27rm//3, and 

0J fc (t) ee cos(wJft) exp(- 7 £t), (34a) 
(*) = IhoySt + sin^f )] exp(- 7 £f ) . (34b) 

The 77-coefficients in Eq. (|33|l are all complex in general 
due to the FDT; except those arising from the Mat- 
subara contribution are real in a canonical bath ensemble; 
see Eqs. (|X6 ]) -(|X9 j) . 

Note that <j%(t) = e _ To* and ffi(t) = ^te-^ 1 , as 
lu% = 0. Also (noting uj' k a = 4o7o + u k ) 

9t4Qk(t) = -7»W - "Ztfk+M ( 35a ) 

d t ^ k+1 {t) = < a ^ k (t) - Tjf^k+iW- (35b) 

The above closed relations will be used in the following 
IGF-COPI construction of hierarchical EOM. 

In principle, the parametrized C a in Eq. (|33p can be 
exact for arbitrary dissipation, if K and M in Eq. (|33[) 
are large enough. The hierarchical EOM to be developed 
in the rest of this section will also be exact; but its size 
grows in a power law. The exact evaluation of complex 
dissipation would rapidly become extremely tedious. We 
will come back to this issue on how to incorporate the 
residue correction, due to the small difference between 
the exact C a (t) (including the zero temperature case) 
and the practically used ones, into the final theory in 
Sec.El 

B. Hierarchical construction: Influence generating 
functionals 

Note that Q a [Eq. (fTTj) ] appears additive with respect 
to the individual components of C' a . The dissipation 
functional TZ [Eq. (|13b|) with Eq. (JTDJ)] will also be addi- 
tive. Moreover, in comparison with the first two terms in 
Eq. (j33j) . the Matsubara term possesses the special prop- 
erties for its pre-exponential factors fj^ being real and its 
time constants j m being dissipation- mode independent. 
As results, the Matsubara contributions to TZ [Eq. (|13bj) ] 



can have the summation over a' performed in the follow- 
ing construction of the hierarchical EOM formalism. 
To proceed, let us denote [</>™(i) = exp(— 7™i)] 

Qf(t;{a}) = f dr^{t-T)Q a ,[a{T)], (36) 

Jto 

and 

Q a m (t; {a}) =J2^ [ dre-^-^QAair)}. (37) 

The corresponding composite Q-functionals that specify 
the dissipation functional [cf. Eq. (|13b[) ] are denoted as 

X£(t; {a}) ee - V2k*QTk) > (38a) 

Y k a (t; {«}) = -zfe+i^fc+i - %Y + i<§2 a fe+ i)>(38b) 

ZS(t; {a}) ee - <*Q'£), (38c) 

and [noting that 77^ is real; see Eq. (|A9|) ] 

Z»(t; {a}) = (38d) 

Included in each of the equation is also the factor of (— i) 
for the sake of bookkeeping; e.g., Z£ amounts to the 
-iQ a in Sec.[mi 

The dissipative functional TZ, by which 

d t T = -1ZT, (39a) 

reads now [cf. Eq. l[Tg])] 

TZ = iJ2 QaZS + * E Q ^ + Y k) + * E Q A 

(39b) 

Apparently, all composite Q-functionals, Eqs. (j3"8]) , are 
influence generating functionals; they are however not 
completed. 

The crucial step in the IGF-COPI construction of hier- 
archical EOM is the time derivatives on these composite 
Q-functionals (cf. Sec. lIIip . Unlike the Z-functionals for 
the Drude and Matsubara components, the time deriva- 
tives of the X- and F-functionals are closed together with 
two additional non-composite Q-functionals [cf. Eq. (|35|) ]. 

X%(t; {a}) ee -i(r% k Q% k+1 - VskQTk+i), (40a) 
Y k a (t; {«}) ee -i« fe+1 QJ fe - V^+iQTk)- (40b) 

The time derivatives of all involving Q-functionals are 
obtained as [cf. Eqs. (35]) and ([21]) ] 

d t X% =A%- ~, a k Xt - L0%Xg, (41a) 
d t X% = J k a X% - 7 fcA£; (41b) 

d t Y k a = B^-^Y k a -^Y k a , (41c) 
d t Y k a =u;' k a Y k a -^Y k a ; (41d) 
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and 



d t ZS = CZ- 7 £ZS, d t Z a m =C a m --f m Z a m . (41c) 



Here 



At 

na 



(42a) 
(42b) 
(42c) 
(42d) 



Not e tha t Eq. ([42cj) is identical to Eq. ([25]). In writing 
Eq. (|42dp , we have used the property that f]^ is real; see 
Eq. JX9]). 

The above six (XKZ)-functionals constitute now a 
complete set of IGF. The general expression for the AIFs 
in the hierarchy are then obtained as [cf. Eq. (|27[) and its 
comment above] 

a,k 

xH(zS) n? >l[(z?nf m }^ (43) 

The index in J- n is specified by the involving nonnegative 
integers, 

n = (n",n",n",n^; j = 0,- • -,2K + 1; m = l,---,M) ■ (44) 

Therefore, the total number of the nonnegative integers 
in the index-set n is 4(K+ l)p + p + Mq, where q denotes 
the number of dissipative modes, and p < q 2 the number 
of nonzero dissipative mode pairs. 



C. Hierarchical equations of motion 

The IGF-COPI approach to the hierarchical EOM can 
now be completed by taking the time derivative on T n 
[Eq. (g3]) with Eqs. ([S2]) and (gT])]. The final results read 
in terms of the auxiliary density operators as [cf. Eqs. (|29|) 
and (J3UJ)] 



„{-} i „{+> 
Pn ~r Pn 



(45) 



The 7-term in Eqs. (|45p arises from the damping-terms 
of Eqs. (|41[) . The resulting damping constant is given by 



7n = E( n "fc + U 2k + n 2k+l + ~ n 2k+ 
a.k 



(46) 



The second term in Eq. (|45"1) stems from the (off- 
diagonal) swap-terms of Eqs. (|41a)l — (|41d|) . It reads 



Pn 



= -J2 UJ k{n2kPn^ 2k +n2k+lPn- 2 



a.k 



Y,^ a ( fl 2kPn- 2k +^ k+1 p nZ2k + 1 ). (47) 



The index-set 
(»*,n*)-(n£ 



n^- differs from n of Eq. (|44|) only by 
- 1, m + 1), while by (n" + 1, nf - 



1) <— (n",n™), at the specified (a,j). 

The third term in Eq. (g5]) stems from the (A, B, C)- 
terms of Eqs. (|41l) . while the last term is from Eq. (|39|) . 
They are the hierarchy-down and hicrarchy-up contribu- 
tions, respectively, and given by 

pt } = E (n% k A%p n - 2k + n a 2k+1 Btp,- 2k+ ) 

a,k 

+ E D + E *Aptz, m > (48) 



and 



Pn 



D,0, 1, ■ 



-Up u 



,2K + 1) 



i E Sa ^n+ . - * E 2a ^n+„ 



(49) 



The index-set ■ ■ or n^) differs from n only by 
changing the specified n™ (fi™ or nj^) to n" ± 1 (fi° — 1 
or ±1). Note that the (n+ -^-auxiliary reduced 
density operators are not generated from Eq. (|4*5j) . since 
Xk and Yfc do not appear in the dissipation functional 1Z 
[Eq. ([39]) ] ; they are rather generated from the EOM for 
the (n+ -^-auxiliary reduced density operators via the 
involving swap {^}-terms there [cf. Eq. (|47|) ]. 

In Eq. (gH), A%, B", C", and denote the reduced 
Liouville-space operator counterparts of Eqs. (|42|) [cf. the 
comments above Eq. (|32p]. 



^0 = 

c°d = 



-i(rtQa'd- v z*dQ a ,), 



(50a) 
(50b) 
(50c) 
(50d) 



a.k 



The reduced Liouville-space operator Q a involved in 
Eq. igSJ) was given by Eq. ©. 

The initial conditions to Eqs. (|45|) are /o n (*o) = 
p(to)5 n o, as inferred from their definitions, and po{i) = 
pit) is the reduced density operator of primary interest. 
Note that when the initial time to — > — oo, the estab- 
lished Hierarchical EOM formalism imposes no approx- 
imation. The initial conditions are however p n (to) = 0, 
where io can be any time before applying the external 
time-dependent fields. The pulse-field induced dynamics 
will then be evaluated via Eq. (|45|) . 

The hierarchical EOM, Eqs. (|I5 |) -([50] l . are exact for a 
general dissipation system that involves the parametrized 
bath correlation functions of Eqs. ((33)) . The residue ef- 
fect due to the small difference between the exact and the 
parametrized C a (t) on the final formalism will be carried 
out together with the hierarchy truncation via the prin- 
ciple of residue correction in the coming section. 

It is noted that with a proper re- arrangement, Eqs. (|45j) 
can be recast in the standard tridiagonal coupling form; 
see Appendix[B] Included there is also a variation of 



7 



the above theory, expressed in terms of the hierarchical 
Green's functions and their related memory kernels and 
continued fraction formalism. 

V. TRUNCATION AND THE PRINCIPLE OF 
RESIDUE CORRECTION 

A. The principle of residue correction 

To complete the EOM formalism, the infinite hierar- 
chy in Eqs. (|45|) - ([50|) should be truncated at a certain 
level, say the (-/Vt run ) th -tier. The simplest way is to set 
all {/?„} of the higher tiers to be zero. The resulting 
Po{t) = p(t) of the primary interest will be exact up to 
the (2A r trun) th -order in the system-bath coupling. Other 
truncation schemes related to different ways to the par- 
tial account for the higher orders effect, have also been 
proposed [U GJ, M, El, El . 

Here, we present the principle of residue correction, 
which itself is formally exact. It is applied not just to 
the truncation, but also to a recovery of the residue ef- 
fect, due to the difference between the exact and the 
parametrized C a (t) of Eq. (f33|) . at all levels of hierarchy. 
The hierarchy truncation that concerns only at the an- 
chor level will be treated in the next subsection. 

The principle of residue correction related to the finite 
difference between the exact and the parametrized ones, 

5C a (t) = C c ™(t)-C a (t), (51) 

arises from the observation that the dissipation func- 
tional is additive [cf. Eq. (fT3b|) with Eqs. (JTUD and (TTTj) ]. 

K cx [t; {a}} = TZ[t; {a}} + 5K[t; {a}]. (52) 

Here, 

8n[t;{a}]=J2Qa[a]5Q a (t;{a}), (53) 

a 

with SQ a the same as Eq. (flUf . but associating with 
SC a (t); i.e. [cf. Eq. (fITa|)] 

6Q a (t;{a})= f dT6C a {t-T)Q a ,[a{T)}. (54) 

The resulting exact influence functional of primary in- 
terest reads now ^ rex [a] = ■F re si[a:].7 r [a!], with dtJ- IC si = 
-STZTrcsi. The AIFs defined in Eq. (gSJ) for the construc- 
tion of the hierarchical EOM should now be replaced by 

3" n = -F r\3~ resi- (55) 

Its time derivative reads 

dtK™ = (^.Frcsi - SUP:™. (56) 

The EOM for the corresponding exact p n is then obtained 
as [cf. Eqs. (|i5 |) -([513 ]) ] 

Pn = - [iC + 5K(t) + 7n ]p n + P r } +P^ } +pi +} - (57) 



The residue correction due to 5C a (t) [Eq. (jSTj) ] is thus 
global; the resulting 51Z(t) [Eq. ([53")) ] modifies the indi- 
vidual hierarchical EOM at all levels. 

Apparently, 6K of Eq. (53]) and K of Eq. (fl3b|) are 
of the same mathematical structure. If their common 
operator-level expression were known and implementable 
readily without approximation, Eq. (|14b|) would be used 
directly, without invoking the hierarchical EOM at all. 
However, it is only possible for special cases, such as the 
pure-dephasing limit (i.e., the case ofj-ff, Q a ] = 0) or the 
driven Brownian oscillator system [16j |. 

The key idea behind Eq. ([57)1 for the general non- 
Markovian dissipation is as followings. The total C™ a is 
partitioned into two parts. One is the parametrized C a 
that carries most of the non-Markovian coupling strength 
and is expressed in the form of Eq. (f3"3"| . Another is 
the residue 6C a that is assumed in the weak interaction 
regime. The strong dissipation due to parametrized C a is 
treated via the hierarchical EOM approach developed in 
Sec. lIVI without approximation. In principle, the residue 
5C a can be zero if K and M for the parametrized C a in 
Eq. (|33p are large enough. However, the size of hierarchi- 
cal EOM grows in a power law as K and M increase; the 
exact evaluation of complex dissipation would rapidly be- 
come extremely tedious if not impossible. Therefore, it is 
a practical trade-off to have a nonzero but weak 5C a , as 
long as its induced global residue 51Z can be accurately 
described with a certain perturbative or nonperturbative 
formulation at the operator level. 

Let us start with the simplest one, the Markovian- 
rcsiduc limit, in which Eq. (1541) reduces to 

SQ a ™ SC a (t)Q a ,[a(t)], (58) 

with 

SC a (t) = f dr5C a {t- T ). (59) 

Jto 

Note that the system variable Q a > in Eq. (|58|) is now rep- 
resented at the ending time t of the path integral at which 
a(t) = a and a'{t) = a' are fixed. As results, Eq. l|5"3"|) 
for the Markovian-residue dissipation can be expressed 
in the operator level as 

5K{t)6 = J2iQa,SC a (t)Q a ,d ~ [5C a {t))*6Q a ,]. (60) 

a 

When t — > co, the above equation reduces to the Matsub- 
ara residue or finite temperature correction proposed by 
Ishizaki and Tanimura [111, ■ Their hierarchical EOM 
is the single- mode Drude dissipation version of Eq. (|57]) . 

The principle of residue correction is right rooted at 
the fact that 51Z [Eq. (|53")) ] is of the same mathemat- 
ical structure as time-local dissipation functional 1Z(t) 
[Eq. (|13bp ]. As results, various well-established approx- 
imation schemes can be exploited for the superoperator 
51Z, as it describes weak residue dissipation. The most 
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celebrated scheme may be the second-order time-local ex- 
pression [11, Bill El, 

5K(t)0 « £ [Q^Q { a ] VP ~ d[6Q^(t)V ] , (61a) 

a 

where 

6QW= [ dTSC a (t-T)e- iC ^Q a ,. (61b) 

J t 

The dissipation-free propagator (assuming also time- 
independent system Hamiltonian for simplicity) is used 
here to connect Q a '[ a ( T )] i n Eq. (|54|) to its value at the 
fixed ending path point of a(t) — a. The above oper- 
ator level expressions are thus obtained [~17^ . In fact, 
Eqs. (fBTjl constitute the unified Bloch-Redfield-Fokker- 
Planck equation [3, [l6[ . 

Nonperturbative approximation schemes can also be 
applied to the evaluation of the global residue-induced 
SlZ(t). These include the non-interacting blip approxima- 
tion or its variations [!, d, [HJ , the self-energy augmented 
methods [35[ , and the use of the exactly solvable driven 
Brownian oscillator model fl6j to mimic the anharmonic 
system of interest. 

B. Hierarchy truncation via the principle of 
residue correction 

Consider now the hierarchy truncation and its related 
issues. Apparently, the collection of anchor {/?n} should 
be properly specified. The hierarchy tier-up p^ +> associ- 
ating with an anchor p^ contains at least one component 
that goes beyond the desired anchoring confinement, and 
thus, is subject to truncation. 

The principle of residue correction is applied here by 
recognizing that the tier-up components can be recast as 
[cf. Eqs. O] 

Pn+ +Pn+ =(*fc+n>N; (62a) 

a, 2k a,2k-\-l 

Pn+ d = z Spn; p N + m =Z^p N . (62b) 

The (2fc) th - and (2k + l) th -components for each dis- 
sipative mode pair a are grouped together for the 
reasons that they arise from the same term of the 
parametrized spectral density function and they carry 
the same strength; cf. Appendix[X] Therefore, they shall 
be treated equally as the truncation is concerned. 

The anchoring indexes can now be specified for the 
individual constituents of the interaction bath correlation 
functions C a (t) of Eq. (|33[) as 

NS, AC; with *; = (),••• ,K; m=l,---,M. 

The closed set of hierarchically coupled EOM will then 
contain those p n , whose individual index-set consists of 
the nonnegative integers that are confined within 

^2k + «2 fe+ i + "2k + nZk+i < N k > (63a) 
<<N£, <<A«. (63b) 



The anchor index-set N in p^ can now be defined as those 
with at least one of the above upper limits being reached. 
The constraint of Eq. (|63a[) is consistent with the way of 
grouping in Eq. (|62a[k Once the upper limit of Eq. (|63aJ) 
is reached, the associated tier- up p w + and are 

"a,2k "a,2k + l 

both subject to the truncation. 

The truncation can now be made based on the formally 
exact relations of Eqs. The involving X%, Y k a , Z£ 

and Zf n there are the reduced Liouville-space oper- 
ators, whose path-integral representation counterparts 
were given by Eqs. ((Ml). Taking Eq. (|3"5aj) for X% for 
an example, its operator-level form reads in contact with 
Eq. (|62a|) as 

p Ht 2k EE X%p n = -l[^ k QUt)PK-ll2k*PK #)]• (64) 

If /°n+ goes beyond the closed hierarchy, an approxi- 
mated expression of Q% k (t) , such as [cf . Eq. (|61b|) and the 
comments there] 

QMQ « I dT<tq k {t ~ r)e- lC ^Q a ,, (65) 

is adopted locally in the rhs of Eq. (|64[) to make the trun- 
cation. The resulting p N + retains the same leading 

a, 2k 

term as the exact one, being of (2A + 2) th -order in the 
specified system-bath coupling strength. Thus, the trun- 
cation with sufficiently large N induces practically no er- 
ror as far as the dynamics of p(t) = po(t) of primary 
interest is concerned. The construction of closed hierar- 
chical EOM with residue correction is now completed. 

C. Discussions and comments 

We re-emphasize here that when the values of K, M 
in Eq. (|33|) and the truncation anchor indexes in Eq. (|63p 
are set to be sufficiently large, the residue effects on the 
primarily interested p are effectively zero. The global 
residue correction and the truncation scheme introduced 
at both the global and the local truncation levels in the 
previous two subsections are made for the purpose of ef- 
ficient evaluation of the reduced dynamics of primary in- 
terest. For example, one may simply terminate the hier- 
archical EOM by setting the aforementioned beyond-the- 
hierarchy p N + = 0. The resulting p of primary inter- 
est will be exact up to the (2A) th order in the specified 
system-bath coupling strength. The improved trunca- 
tion as Eq. (B3]) will lead to p exact up to the (2A + 2) th 
order, rather than the (2A) th order. The above two trun- 
cation schemes, which represent two different resumma- 
tions for partially incorporating the higher-order effects 
on the reduced dynamics of primary interest, shall be of 
no practical difference when the convergence is reached. 

Note that various commonly used forms of the Bloch- 
Redfield theory and Fokker-Planck equations can be con- 
sidered as the globally weak (residue) dissipation with- 
out invoking the hierarchical EOM at all. They can also 
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be recovered if the second-order truncation scheme of 
Eq. (|6"5j) is applied to the primary tier of N = 0. On 
the other hand, if all second-tier auxiliary reduced den- 
sity operators are set to be zero, the present hierarchi- 
cal EOM formalism reduces to the second-order memory- 
kernel quantum dissipation theory. 

The global residue dissipation 81Z(t) introduced in the 
final formalism [Eq. ([57)) ] is also for the purpose of ef- 
ficiency. Note the number of the n th -tier auxiliary re- 
duced density operators {p n } is , where P = 
4(K + l)p + p + Mq is the number of the nonnegative 
integers in the index n; see comments after Eq. (14"4"| and 
in AppcndixlBl before Eq. (|B 1 1) . The size of the hierar- 
chical EOM increases in a power law as the values of 
K and M for the parametrized C a of Eq. (|33[) . The 
residue correction is introduced to reduce the required 
values of K and M, as long as the induced residue weak 
dissipation 51Z(t) can be accurately evaluated. In this 
sense, the final residue-corrected hierarchical EOM for- 
malism remains practically exact. It is worthy to point 
out that the final formalism is even capable of treating 
the zero-temperature dissipation. At T = 0, the Matsub- 
ara frequencies are all vanished and the exact C° xa (i) is 
given by Eq. (|A10[) . However, one can set a certain finite 
low temperature for the parametrized C a , and evaluate 
residue (5C a -corrected hierarchical EOM dynamics, fol- 
lowed by the convergence test with a lower value of the 
parametrization temperature. 

VI. SUMMARY 

In summary, we have constructed the residue-corrected 
hierarchical EOM formalism [Eq. dSTJ) with Eqs. (|4l)j)- 
[5D])]. The construction consists of two main steps. The 
first is the IGF-COPI approach to the EOM formalism 
fSec. HVI) . for the parametrized interaction bath correla- 
tion functions preserving the fluctuation-dissipation the- 
orem. This step is itself exact, as the FDT-preserved 
parametrization involved can in principle represent ar- 
bitrary interaction bath correlation functions. The sec- 
ond step is the residue correction, concerning about the 
practical applicability of the present theory to a broad 
range of systems. The principle of residue correction 
(Sec.|V| is itself exact, rooted from the formal relation 
between the dissipation functional and the time-local dis- 
sipation superoperator; see Sec. Ill CI The application of 
this principle to construct the global (5C-induced) or the 
local (truncation-induced) residue-correction invokes in- 
evitably a certain approximation scheme. However, it 
can be made in a sufficiently accurate manner, as far as 
the primarily interested p(t) is concerned (cf. Sec. IV C[) . 
As results, the residue-corrected hierarchical EOM for- 
malism could be practically used in the study of general 
quantum dissipation systems interacting with arbitrary 
bath, arbitrary time-dependent external fields, and at 
any temperature, including T — 0; see the last remark 
stated in Sec. lVCl 



The hierarchical EOM formalism may be relatively 
tractable, in comparison with the direct evaluation of the 
path-integral formalism jlg| . The memory effect in the 
quasiadiabatic propagator path integral method is de- 
scribed in terms of the nonlocality of the influence func- 
tional [18j |. while it in the present differential formalism 
is resolved via a set of linearly coupled time-local aux- 
iliary operators. The exponential-like series expansion 
of the parametrized C a (t) [Eq. (|33)) ] can be considered 
as the separation of the time scales, and the resulting 
p n is associated with the decay constant 7 n . The dif- 
ferent truncation anchor indexes [Eq. (|63p ] can therefore 
be identified to reduce the required number of equations, 
which otherwise grows exponentially if all C a -composites 
are treated equally. The weak residue correction at the 
global level provides an additional freedom to improve 
the numerical efficiency of the present EOM theory. 

Numerous existing quantum dissipation theories can 
be recovered from the present formalism. These in- 
clude the Tanimura's hierarchical EOM for single-mode 
Drude dissipation [2l|, [22j and the unified Bloch-Redfield- 
Fokker-Planck formulation [3, EH; see comments af- 
ter Eqs. (jnOl) and the second paragraph of Sec. lVCl re- 
spectively. The equivalent hierarchical and continued- 
fraction Green's functions and memory kernels expres- 
sions are also presented; see Appendix[B] The applica- 
tion of the continued-fraction Green's function formalism 
to the two-state electron transfer system, with a single- 
mode Drude-Debye dissipation at the high-temperature 
limit, has been carried out recently, resulting in an an- 
alytical expression for the nonperturbative rate process 

IE 13. 

Multi-mode dissipation is physically important. In the 
weak (second-order) dissipation regime, the effects of dif- 
ferent system-bath coupling modes are additive. This 
simple property is no longer true in the strong dissipation 
regime, and it has been shown that the strong dissipation 
should include also the cases of long memory system- 
bath interactions [l7l |. Cooperative dissipation, similar 
like the co-tunneling in quantum transport [34j , could be 
a common phenomenon in reality. The present work is 
carried out in the bosonic canonical bath case. The ex- 
tension to the grand canonical bath ensemble cases, in- 
cluding both fermion and boson statistics, will be treated 
elsewhere |33| . 
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APPENDIX A: PARAMETRIZATION OF BATH 
CORRELATION FUNCTIONS 

The fluctuation-dissipation theorem (FDT) relates the 
correlation functions C aa ' if) and the spectral density 
functions J aa >(uj) by 



C aa> {t) 



duj- 



1 



(Al) 



The spectral density functions J aa i (u>) in a canonical en- 
semble satisfy in general the symmetry relations 



Jaa'(w) = -Ja'a(-W) = J^ a (w). 



(A2) 



The FDT leads also to J aa '{0) — 0, the spectrum pos- 
itivity J ao (w > 0) > 0, and the Schwarz- inequality 

Jaa(oj)Jbb(u) > \Jab{u)\ 2 - 

We adopt the following form of the extended Meier- 
Tannor spectral-density parametrization scheme [l|| [38| , 
in which (setting = 0) 



A' 



J(u) - ,f- C fc a 7 fc ^ + <> 2 fA3) 

Ja[ - UJ ) - ,.,2 , („a\2 + |,.,2 _ (,.,a , ,-^212 ' 



All involving parameters are real; wj^g, 7™ and 7° posi- 
tive as well. The symmetry relations in Eq. (|A2p require 
also (noting that C£ a = 0) 

/ „,ba /-ba , ,ba „ba yba 7ba\ („,ab /-ab , ,ab „,afc /-ab 7-ab\ 
17d iCn > w fe >7fc ,Cfc ,Cfc J = (,7d >Cd ,7 & ,C fe ,-Cfc J- 

The FDT [Eq. ([AT]) ] leads to the interaction bath corre- 
lation functions of [cf. Eq. (B4) of Ref. EH 



2K+1 

c a (t>o) = < e -^+ $>M(*)+ 
3=0 



rn — 1 



-a g-7m« 



(A4) 



with (noting u>q = 0) 

^ fe (i) = cosKt)exp(- 7 ^), 



(A5a) 



4>2k+i (*) = [4o7o * + sin(w|ft)] exp(-7f *)« (A5b) 



The first term in the rhs of Eq. (|A4[) arises from the 
pole of the Drude term in Eq. (|A3[) . The involving coef- 
ficient is given by 



Sd 



ctan(/3 7 »/2)-i] 



(A6) 



The second term in the rhs of Eq. (|A4[) arises from 
the poles of the second term in Eq. (|A3[) . The involv- 
ing functions (j>j(t), given by Eqs. (fM)) . are chosen to be 
real in this work, rather than complex exponential func- 
tions used in Ref. |lg. Noting that </>g(i) = e _7ot and 
= j£t(j)%(t), as inferred from = 0. The in- 
volving coefficients {rj°'} are complex. Those with even 
indexes are (including k = at which ujq = 0), 



4< 7 «[cosh(/30-cos(/3 7 -)] 



,a 

4 7 £' 



(A7a) 



and those with odd indexes are 

.a (Co a + Co) 



ctau087ff/2)-i] 



47o" 

and (2fc + 1 = 3,5, • • ■) 

C>£ sinh(/3^) + (Q + C£) 7 £ sm(P^) 



(A7b) 



V2k+1 



4w»7»[cosh(/3w°) - cos(/3 7 ^)] 
4w? 



(A7c) 



The last term in the rhs of Eq. (|A4p arises from the 
Matsubara frequencies, 



7 m = 27rm//3; m = 1,2, ••• , 



(A8) 



which constitute poles of the (1 — e 13 z ) 1 factor in 
Eq. (|Aip . The involving coefficients are 



-i(2/P)J a (-i% 



1,2, 



(A9) 



Note that in general, while J a = J aa ' (w) can be complex 
when a/a', its analytical continuation to J a (iuj) is pure 
imaginary in a canonical ensemble system, as inferred 
from Eq. |[A2]) . As results, of Eq. (fAUj) are all real. 

In principle, Eqs. (|A4|) — (|A9|) can be exact at an ar- 
bitrary finite temperature, if K — > 00. In the hierar- 
chical construction presented in Sec. HVl K is finite and 
the Matsubara terms m = 1, • • • , M is also finite. The 
residue correction due to the small difference between the 
exact C a (t) and the parametrized ones with finite terms 
is considered in Sec.fV] By doing that, we can even treat 
the dissipation at T = 0, where the FDT of Eq. (|ATj) 
assumes 



C a (t) 



1 f 00 

- / duje- iuJt J a (uj); for T = 0, (A10) 
I" Jo 



while the parametrized ones used in Sec. lIVI are at a cer- 
tain finite low temperature. 



APPENDIX B: RECURSIVE GREEN'S 
FUNCTIONS, MEMORY KERNELS, AND 
CONTINUED FRACTION FORMALISM 

1. The EOM in tridiagonal coupling matrix form 

For the case of a single-mode Drude-Debye dissipa- 
tion at high temperature, the EOM assumes the stan- 
dard tridiagonal coupling matrix form; i.e., Eqs. (|3ip 
in the special case of single mode. This property has 
been used to construct the continued fraction formal- 
ism [lj| [2(J Hi], [22j]. The involving Green's functions 
are also identified, together with an interesting applica- 
tion to the establishment of an analytical expression for 
electron transfer rate processes (36j . 



11 



Here, we would like to extend our previous continued 
fraction-Green's function formalism to the present com- 
plex dissipation systems. To that end, we shall first re- 
cast Eqs. (f4"5j) in the standard tridiagonal coupling ma- 
trix form. This is done by considering the order of p n 
depending on the hierarchy generating functionals. 



Kfe + ™2fc+l + ™2fc + «2fe+l) + < 



Collect now all those n th -tier auxiliary density operators; 
i.e., p n that have n n = n, and arrange them into a vector 
p n = {p n ; n„ = n}. The size of this vector is ^rpE^r i if 
there are no truncations involved; see Sec. IV Bl Here P 
is the number of nonnegative integers in the index-set n; 
see the comments after Eq. (|4"4")l . 

The hierarchical EOM [Eqs. ([IS])] can then be re- 
arranged in the standard tridiagonal matrix form as 



dtP n = ~A„p„ - iA n p n _ 1 - iB n p n 



(Bl) 



Note that the swap ones p { n ^ } are now a part of p n , and 
the A„ matrix elements are all numbers, except the Li- 
ouvillian £, which can also be time dependent in the 
presence of external pulsed fields for example. All p„~ } 
are arranged into p n _i with the Liouville-space opera- 
tors defined in Eqs. (fST))) involved in the elements of A n . 
Those pn +} are now in p n+1 and the elements in B n are 
according to Eq. ([49]) . 

Note that p (t) = p{t) is the reduced density opera- 
tor of primary interest, and Ao = iC. The collection of 
the n th -tier auxiliary density operators have the leading 
contributions of the (2n) th -order in the overall system- 
bath couplings to the p(t) of primary interest. The initial 
conditions to Eqs. ()B1|) are p n (to) = S n op(tQ), as inferred 
from the definition. 



2. Green's functions versus memory kernels 

In terms of the propagators, by which 

p n (t)=U n (t,t )p(t ), (B2) 
Eqs. (|B1[) read [noting that U n (t ,t ) = 5 n o] 

d t U n = -A n U n - iA n U n -i - iB n U n+1 . (B3) 



Introduce now the hierarchical set of Green's functions 
in the reduced system Liouville space as follows. 

U Q (t,t ) = g {t,to), (B4a) 
U n (t,t ) = -if d,Tg n (t,T)A n U n - 1 (T,t Q ), n > 1. 

•I to 

(B4b) 

The initial conditions are Q n (ta,ta) — E Substituting 
Eqs. (HU into Eq. (|B3j) leads then to 



d t &n(t,t ) = -A n g n (t,t ) - / dTll n (t,T)g n (T,t Q ), 

J to 

(B5a) 

with the involving memory kernel of 

n n (t,r) = B n g n+1 (t,r)A n+1 . (B5b) 

In particular, IIo(t, r) = n(t, r) is the primary mem- 
ory kernel, by which (noting that Aq — iC) 



p(t) = -iC(t)p(t) - / dTlL(t,T)p(r). (B6) 

J t 

For the time-independent system Hamiltonian, 
On(t,r) = g n (t- r) and n„(t,r) = II n (t - r). The 
above Green's functions (or memory kernels) can be 
resolved in the Laplace frequenc y domain via the 
continued fraction expression of [36l l37l| 



g n (s) 



s + A n + B n g n+1 (s)A n+1 



(B7) 



Unlike the auxiliary propagators U n that couple with 
both U n +\ and U n -x, the Green's functions g n cou- 
ple only with Gn+i- Thus, the evaluation of the re- 
duced dynamics of primary interest could be more effi- 
ciently carried out in terms of the auxiliary Green's func- 
tions, especially when there is no time-dependent exter- 
nal fields. The continued fraction formalism, combined 
with the Dyson equation technique, has been recently ap- 
plied to the simple electron transfer, a spin-boson system, 
in Drude-Debye solvents [36|, H3] • 
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